Ensayo de Segmentación de imágenes RGB, basado en el ejercicio publicado por Ivan Lizarazo en formato de Notebook.

library ("rasterVis") |> suppressPackageStartupMessages()
library ('OpenImageR') |> suppressPackageStartupMessages()
library ('RImagePalette') |> suppressPackageStartupMessages()
library("raster") |> suppressPackageStartupMessages()
source(file = "indices_vegetacion_RGB.R")

##1. Leer la imagen para analizar

La imagen ha sido descargada del sitio web del Observatorio de la Tierra de la NASA . Ilustra el diverso paisaje agrícola en la parte occidental del estado de Minas Gerais en Brasil.

Se hace la asignación de la ruta para generalizar

# asignar la imagen con su ruta 
# imagen <- "~/Pictures/balandra.jpg"
# imagen <- "~/giserias/conchalito-vuelo2/conchal_10Q.tif"
imagen <- "~/giserias/conchalito-vuelo2/RPubs_crops.jpg"

Se lee la imagen a analizar:

(crops <- stack(imagen)) # es un RasterStack con estructura para índices
## class      : RasterStack 
## dimensions : 960, 1440, 1382400, 3  (nrow, ncol, ncell, nlayers)
## resolution : 1, 1  (x, y)
## extent     : 0, 1440, 0, 960  (xmin, xmax, ymin, ymax)
## crs        : NA 
## names      : RPubs_crops_1, RPubs_crops_2, RPubs_crops_3

Vemos la imagen cargada:

plotRGB(crops, 1, 2, 3) # visualiza un RasterStack
## Warning: [rast] unknown extent

## Warning: [rast] unknown extent

## Warning: [rast] unknown extent

True Color Vegetation Indices

Thanks to Matthias Forkel we can use a nice palette for vegetation indices:

colores <- function(n) {
    .fun <- colorRampPalette(c("chocolate4", "orange", "yellow", "grey", "green", "green3", "darkgreen"))
    col <- .fun(n)
    return(col)
}

Let’s calculate the EG index:

# For the crops image red = 1, green=2, blue = 1.
(egvi <- EGVI(crops, 1, 2,3)) # los índices regresan un RasterLayer 
## Warning: [rast] unknown extent

## Warning: [rast] unknown extent

## Warning: [rast] unknown extent

## Warning: [rast] unknown extent

## Warning: [rast] unknown extent

## Warning: [rast] unknown extent
## class      : RasterLayer 
## dimensions : 960, 1440, 1382400  (nrow, ncol, ncell)
## resolution : 1, 1  (x, y)
## extent     : 0, 1440, 0, 960  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : layer 
## values     : -0.3543046, 1.33945  (min, max)

Visualizamos:

cortes <- seq(-0.35, 1, by=0.1)
cols <- colores(length(cortes)-1)
# plot() despliega un RasterLayer
plot(egvi, col = cols, breaks=cortes, main = "EG Vegetation Index")

Calculamos el índice NGRDI:

# For the crops image red = 1, green=2, blue = 1.
(ngvi <- NGRDI(crops, 1, 2,3))
## Warning: [rast] unknown extent

## Warning: [rast] unknown extent

## Warning: [rast] unknown extent

## Warning: [rast] unknown extent
## class      : RasterLayer 
## dimensions : 960, 1440, 1382400  (nrow, ncol, ncell)
## resolution : 1, 1  (x, y)
## extent     : 0, 1440, 0, 960  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : layer 
## values     : -0.3809524, 1  (min, max)

Visualizamos:

cortes <- seq(-0.35, 1.0, by=0.1)
cols <- colores(length(cortes)-1)
##plot(ndvitrend, 2, col=cols, breaks=classbreaks)
plot(ngvi, col = cols, breaks=cortes, main = "NGRDI Vegetation Index")

2. Image segmentation

Leamos la imagen usando el paquete OpenImageR:

#RGB.Color <- RGB.Color[, , 1:3] # si lee cuatro capas en lugar de tres
RGB.Color <- readImage(imagen) # Large array
# Nota: las funciones readJPGE, readPNG leen a matriz la imagen, tal y como lo hace
# {OpenImageR}::readImage() lee a matriz o array los formatos  .png, .jpeg, .jpg, .tif .tiff

Usemos ahora SLIC para la segmentación de imágenes. Para obtener una explicación detallada del código, consultar la viñeta del paquete.

# superpixel recibe un Large array o imagen en forma de matriz con tres capas de fondo
Region.slic = superpixels(input_image = RGB.Color, method = "slic", superpixel = 80,
   compactness = 30, return_slic_data = TRUE,
   return_labels = TRUE, write_slic = "",
   verbose = FALSE)
## Warning in interface_superpixels(input_image, method, superpixel, compactness,
## : The input data has values between 0.000000 and 1.000000. The image-data will
## be multiplied by the value: 255!
imageShow(Region.slic$slic_data)

3. Data analysis

str(RGB.Color)
##  num [1:960, 1:1440, 1:3] 0.886 0.855 0.843 0.82 0.749 ...

Let’s get a matrix with NGRDI values:

# ngvi
# egvi
# egvi@data@values
# ngvi@data@values
ngvi.mat <- matrix(ngvi@data@values,
            nrow = ngvi@nrows,
            ncol = ngvi@ncols, byrow = TRUE)
# min(ngvi.mat)
# max(ngvi.mat)

La función imageShow() trabaja con datos que están en el rango de ocho bits 0 – 255 o en el rango [0,1] (es decir, el rango de x entre 0 y 1 inclusive). Sin embargo, no funciona con valores NGRVI si estos valores son negativos. Por lo tanto, escalaremos los valores de NGRDI a [0,1].

#rescale between 0 and 1; ver la funcion {OpenImageR}::NormalizeObject() que se usa mas abajo
normalize <- function(x) {
  min <- min(x)
  max <- max(x)
  return(1* (x - min) / (max - min))
}
ngvi.mat_norm <- normalize(ngvi.mat)
imageShow(ngvi.mat_norm)

ngvi.data <- RGB.Color
# ngvi.data[,,1] <- ngvi.mat_norm
# ngvi.data[,,2] <- ngvi.mat_norm
# ngvi.data[,,3] <- ngvi.mat_norm

En la siguiente sección describiremos cómo crear una segmentación de imágenes de los datos del NGRDI y cómo utilizar el análisis de conglomerados para crear una clasificación del terreno

4. Image classification

Aquí se muestra una aplicación de la función superpixels() para la clasificación de la cobertura terrestre de los datos NGRDI.

ngvi.80 = superpixels(input_image = ngvi.data,
   method = "slic", superpixel = 80,
   compactness = 30, return_slic_data = TRUE,
   return_labels = TRUE, write_slic = "",
   verbose = TRUE)
## Warning in interface_superpixels(input_image, method, superpixel, compactness,
## : The input data has values between 0.000000 and 1.000000. The image-data will
## be multiplied by the value: 255!
## The input image has the following dimensions: 960 1440 3
## The 'slic' method will be utilized!
## The output image has the following dimensions: 960 1440 3
imageShow(ngvi.80$slic_data)

Aquí la estructura del objeto ngvi.80.

str(ngvi.80)
## List of 2
##  $ slic_data: num [1:960, 1:1440, 1:3] 226 218 215 209 191 169 151 136 130 127 ...
##  $ labels   : num [1:960, 1:1440] 0 0 0 0 0 0 0 0 0 0 ...

Es una lista con dos elementos, slic_data, una matriz tridimensional de datos de color de píxeles (no normalizados) y etiquetas, una matriz cuyos elementos corresponden a los píxeles de la imagen. El nombre del segundo elemento sugiere que puede contener valores que identifican el superpíxel al que pertenece cada píxel. Veamos si esto es cierto.

sort(unique(as.vector(ngvi.80$labels)))
##  [1]  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
## [26] 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
## [51] 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73

Hay 72 etiquetas únicas. Aunque la llamada a superpixels() especificó 80 superpíxeles, la función generó 72. Podemos ver qué píxeles tienen el valor de etiqueta 0 estableciendo los valores de todos los demás píxeles en [255, 255, 255], lo que se trazará como blanco.

R0 <- ngvi.80
for (i in 1:nrow(R0$label))
    for (j in 1:ncol(R0$label))
       if (R0$label[i,j] != 0)
          R0$slic_data[i,j,] <- c(255,255,255)
imageShow(R0$slic_data)

La operación aísla el superpíxel en la esquina superior derecha de la imagen, junto con la porción correspondiente del límite. Podemos utilizar fácilmente este enfoque para determinar qué valor de ngvi.80$label corresponde a qué superpíxel. Ahora tratemos con el límite. Una pequeña exploración del objeto ngvi.80 sugiere que los píxeles en el límite tienen los tres componentes iguales a cero. Aislamos y trazamos todos esos píxeles coloreando todos los demás píxeles de blanco.

Bdry <- ngvi.80
  for (i in 1:nrow(Bdry$label))
   for (j in 1:ncol(Bdry$label))
     if (!(Bdry$slic_data[i,j,1] == 0 &
     Bdry$slic_data[i,j,2] == 0 &
       Bdry$slic_data[i,j,3] == 0))
       Bdry$slic_data[i,j,] <- c(255,255,255)
Bdry.norm <- NormalizeObject(Bdry$slic_data)
imageShow(Bdry$slic_data)

La figura muestra que efectivamente se han identificado los píxeles límite. Hay que tener en cuenta que la función imageShow() muestra estos píxeles en blanco con un borde negro, en lugar de negro puro.

Habiendo realizado un análisis preliminar, podemos organizar nuestro proceso de segmentación en dos pasos.

El primer paso será reemplazar cada uno de los superpíxeles generados por la función superpixels() de OpenImageR por uno en el que cada píxel tenga el mismo valor, correspondiente a una medida de tendencia central (p. ej., la media, mediana o moda) del superpíxel original.

El segundo paso será utilizar el procedimiento de agrupamiento no supervisado de K-medias para organizar los superpíxeles del paso 1 en un conjunto de grupos y darle a cada grupo un valor correspondiente a una medida de tendencia central del grupo.

Utilizaremos la función make.segments() para realizar la segmentación. El primer argumento de make.segments() es el objeto de superpíxeles y el segundo es la medida funcional de tendencia central. Aunque en este caso cada uno de los tres colores del objeto ngvi.80 tiene los mismos valores, esto puede no ser cierto para todas las aplicaciones, por lo que la función analiza cada color por separado.

La función, escrita por Richard Plant, es la siguiente:

# Identify a measure of central tendency of each superpixel
make.segments <- function(x, ftn){
# The argument ftn is any functional measure of central tendency
   z <- x
# For each identified superpixel, compute measure of central tendency
   for (k in unique(as.vector(x$labels))){
# Identify members of the superpixel having the given label
      in.super <- matrix(0, nrow(x$label), ncol(x$label))
      for (i in 1:nrow(x$label))
         for (j in 1:ncol(x$label))
            if (x$label[i,j] == k)
               in.super[i,j] <- 1
#Identify the boundary cells as having all values 0
      on.bound <- matrix(0, nrow(x$label), ncol(x$label))
      for (i in 1:nrow(x$label))
         for (j in 1:ncol(x$label))
            if (in.super[i,j] == 1){
               if (x$slic_data[i,j,1] == 0 & x$slic_data[i,j,2] == 0 
                  & x$slic_data[i,j,3] == 0)
                     on.bound[i,j] <- 1
         }
#Identify the superpixel cells not on the boundary
      sup.data <- matrix(0, nrow(x$label), ncol(x$label))
         for (i in 1:nrow(x$label))
         for (j in 1:ncol(x$label))
            if (in.super[i,j] == 1 & on.bound[i,j] == 0)
               sup.data[i,j] <- 1
# Compute the measure of central tendency of the cells in R, G, B
      for (n in 1:3){
# Create a matrix M of the same size as the matrix of superpixel values
         M <- matrix(0, dim(x$slic_data)[1], dim(x$slic_data)[2]) 
         for (i in 1:nrow(x$label))
            for (j in 1:ncol(x$label))
# Assign to M the values in the superpixel
               if (sup.data[i,j] == 1) M[i,j] <- x$slic_data[i,j,n]
          if (length(M[which(M > 0 & M < 255)]) > 0)
# Compute the measure of central tendency
            ftn.n <- round(ftn(M[which(M > 0 & M < 255)]), 0)
         else
            ftn.n <- 0
         for (i in 1:nrow(x$label))
            for (j in 1:ncol(x$label))
               if (in.super[i,j] == 1) z$slic_data[i,j,n] <- ftn.n
           }
      }
   return(z)
   }

Aquí está la aplicación de esa función al objeto ngvi.80 con el segundo argumento configurado como mean.

Podemos ver el resultado.

ngvi.means <- make.segments(ngvi.80, mean)
imageShow(ngvi.means$slic_data)

El siguiente paso es desarrollar grupos que representen tipos de cobertura terrestre identificables. En un proyecto real, el procedimiento sería recopilar un conjunto de datos reales del sitio, pero esa opción no está disponible para nosotros. En su lugar, trabajaremos con la reproducción en color real de la foto del ejemplo.

La cobertura del suelo se puede subdividir utilizando K-medias en cualquier número de tipos. Tenga en cuenta que no se ha interpretado aquí la relación entre los grupos de K-medias y las clases de cobertura del suelo existentes en el sitio del ejemplo. Sin embargo, en una aplicación real, esta tarea es imprescindible.

set.seed(123)
# OJO: en la función kmeans se establece el número de clases subsecuentes
nclass <- 10
ngvi.clus <-  kmeans(as.vector(ngvi.means$slic_data[,,1]), nclass)
vege.class <- matrix(ngvi.clus$cluster,
                     nrow = ngvi@nrows, 
                     ncol = ngvi@ncols, byrow = FALSE)
class.ras <- raster(vege.class, 
                    xmn = crops@extent@xmin,
                    xmx = crops@extent@xmax, 
                    ymn = crops@extent@ymin, 
                    ymx = crops@extent@ymax, 
                    crs = crs(crops))

A continuación podemos usar la función {raster} ratify() para asignar niveles de factores descriptivos a los grupos o clases definidas.

class.ras <- ratify(class.ras)
rat.class <- levels(class.ras)[[1]]
# eventualmente, se asignan etiquetas para cada una de las clases del resultado;
# en el caso del ejemplo, se aplica:
# rat.class$landcover <- c("LC 10", "LC 09", "LC 08", "LC 07", "LC 06", "LC 05", "LC 04", "LC 03", "LC 02", "LC 01")
rat.class$landcover <- paste("Clase", nclass:1)
#Crear una palette con el número de colores correspondientes al resultado en kmeans
# se lee de nuevo la imagen original como raster comun en Large array x, y, 1:3
# ncrops <-  jpeg::readJPEG("~/Pictures/balandra.jpg")
ncrops <- readImage(imagen) # para no invocar otra librería
crops.pal <- image_palette(ncrops, n=nclass)
# Nota: las funciones readJPGE, readPNG leen a matriz la imagen, tal y como lo hace
# {OpenImageR}::readImage() lee a matriz o array los formatos  .png, .jpeg, .jpg, .tif .tiff
levels(class.ras) <- rat.class
levelplot(class.ras, margin=FALSE,
   col.regions= crops.pal[1:nclass],
   # se pueden reasignar los colores así como agrupar asignando el mismo a diferentes clases
   # col.regions= crops.pal[c(5,7,9,4,8,10,2,4,1,6)],
    main = "Land Cover Types in crops")